###########################
# This file executes all the procedures to analyze the code and produce tables 
# and figures for "Inaccurate Forecasting of a Randomized Controlled Trial"
# authored by Mats Ahrenshop, Miriam A. Golden, Saad Gulzar, and Luke Sonnet 
# Code written by Mats Ahrenshop
# Reviewed and edited by MG
########################

#---- 0. HOUSEKEEPING ----

library(sandwich)
library(bookdown)
library(tidyverse)
library(lubridate)
library(readxl)
library(estimatr)
library(stargazer)
library(kableExtra)
library(texreg)
library(gridExtra)
library(RColorBrewer)
library(ggpubr)
library(gt) # go get this run `remotes::install_github("rstudio/gt")`
# You might need to `install.packages("remotes")` first
library(modelsummary)
library(nnet)
library(interactions)
library(huxtable)
library(emmeans)
library(tabulator)
library(gmodels)
library(here)
library(readxl)
library(ggpubr)
library(data.table)

all <- read_csv("../dataverse/forecasting_clean.csv")
ivr_mus <- read_csv("../dataverse/ivr_mus.csv")
tab_names <- excel_sheets(path = "../dataverse/rct_counts.xlsx")
list_all <- lapply(tab_names, function(x) read_excel("../dataverse/rct_counts.xlsx", sheet = x))


#---- 1. MAIN PAPER ANALYSIS ----

#### table 1 factorial experiment ####

treatments <- all %>%
  group_by(pilot_first, compliance_first) %>%
  summarize(n = n()) %>%
  filter(!is.na(pilot_first))

treatments

#### table 3 p-values ####

forecast_mus <- all |>
  select(call_eval_mpa, call_eval_gov, call_account,
         responsive_eval_mpa, responsive_eval_gov, responsive_account,
         ans_phone, ans_quest) |>
  summarize_all(mean, na.rm = TRUE)

forecast_mus <- as_vector(forecast_mus)

ivr_mus <- ivr_mus |> add_column(forecast_mus)

ivr_mus

## get p-values for table
# pilot
t.test(all$call_eval_mpa, mu = ivr_mus$pilot_mu[1])
t.test(all$call_eval_gov, mu = ivr_mus$pilot_mu[2])
t.test(all$call_account, mu = ivr_mus$pilot_mu[3])
t.test(all$responsive_eval_mpa, mu = ivr_mus$pilot_mu[4])
t.test(all$responsive_eval_gov, mu = ivr_mus$pilot_mu[5])
t.test(all$responsive_account, mu = ivr_mus$pilot_mu[6])

# intervention
t.test(all$call_eval_mpa, mu = ivr_mus$intervention_mu[1])
t.test(all$call_eval_gov, mu = ivr_mus$intervention_mu[2])
t.test(all$call_account, mu = ivr_mus$intervention_mu[3])
t.test(all$responsive_eval_mpa, mu = ivr_mus$intervention_mu[4])
t.test(all$responsive_eval_gov, mu = ivr_mus$intervention_mu[5])
t.test(all$responsive_account, mu = ivr_mus$intervention_mu[6])
t.test(all$ans_phone, mu = 73.1)
t.test(all$ans_quest, mu = 23.8)

#### table 4 effect of treatments on ITT forecasts ####

# CALL AVG
call_avg_p <- lm_robust(call_avg ~ pilot_first, data = all)
call_avg_c <- lm_robust(call_avg ~ compliance_first, data = all)
call_avg_pc <- lm_robust(call_avg ~ pilot_first*compliance_first, data = all)

## for appendix:
# CALL EVAL MPA
call_mpa_p <- lm_robust(call_eval_mpa ~ pilot_first, data = all)
call_mpa_c <- lm_robust(call_eval_mpa ~ compliance_first, data = all)
call_mpa_pc <- lm_robust(call_eval_mpa ~ pilot_first*compliance_first, data = all)

# CALL EVAL GOV
call_gov_p <- lm_robust(call_eval_gov ~ pilot_first, data = all)
call_gov_c <- lm_robust(call_eval_gov ~ compliance_first, data = all)
call_gov_pc <- lm_robust(call_eval_gov ~ pilot_first*compliance_first, data = all)

# CALL EVAL ACC
call_acc_p <- lm_robust(call_account ~ pilot_first, data = all)
call_acc_c <- lm_robust(call_account ~ compliance_first, data = all)
call_acc_pc <- lm_robust(call_account ~ pilot_first*compliance_first, data = all)

# RESP AVG.
resp_avg_p <- lm_robust(resp_avg ~ pilot_first, data = all)
resp_avg_c <- lm_robust(resp_avg ~ compliance_first, data = all)
resp_avg_pc <- lm_robust(resp_avg ~ pilot_first*compliance_first, data = all)

## for appendix:
# RESP EVAL MPA
resp_mpa_p <- lm_robust(responsive_eval_mpa ~ pilot_first, data = all)
resp_mpa_c <- lm_robust(responsive_eval_mpa ~ compliance_first, data = all)
resp_mpa_pc <- lm_robust(responsive_eval_mpa ~ pilot_first*compliance_first, data = all)

# RESP EVAL GOV
resp_gov_p <- lm_robust(responsive_eval_gov ~ pilot_first, data = all)
resp_gov_c <- lm_robust(responsive_eval_gov ~ compliance_first, data = all)
resp_gov_pc <- lm_robust(responsive_eval_gov ~ pilot_first*compliance_first, data = all)

# RESP ACC.
resp_acc_p <- lm_robust(responsive_account ~ pilot_first, data = all)
resp_acc_c <- lm_robust(responsive_account ~ compliance_first, data = all)
resp_acc_pc <- lm_robust(responsive_account ~ pilot_first*compliance_first, data = all)

texreg(l = list(call_avg_p, call_avg_c, call_avg_pc,
          resp_avg_p, resp_avg_c, resp_avg_pc),
       include.ci = FALSE,
       digits = 3,
       booktabs = TRUE,
       label="tab:main-results",
       stars = c(0.01, 0.05, 0.1),
       caption = "Experimental effects for ITT forecasts (average across indices)",
       caption.above = TRUE,
       custom.header = list("Forecast call on ITT" = 1:3,
                            "Forecast responsive on ITT" = 4:6),
       custom.model.names = c("(1)", "(2)", "(3)", "(4)", "(5)", "(6)"),
       custom.coef.map = list("(Intercept)" = "Constant",
                              "pilot_first" = "Pilot prime",
                              "compliance_first" = "Compliance prime",
                              "pilot_first:compliance_first" = "Interaction P x C"),
       reorder.coef = c(2, 3, 4, 1),
       include.rsquared = F, include.adjrs = F, include.rmse = F, include.nobs = T,
       custom.note = "\n\\item %stars. OLS estimates with HC2 standard errors. \\textit{Forecast call on ITT} refers to the comparison of the treatment arm in which households received a call with a message and a question from their MPA versus the control group that received no call. \\textit{Forecast responsive on ITT} refers to a comparison of the treatment arm in which households received a call with  a message, a question, and a responsive follow-up call from their MPA (mentioning the action taken by the MPA in response to feedback) versus the arm  that received a generic follow-up call (thanking citizens for the input).",
       float.pos = "h",
       #dcolumn = TRUE,
       threeparttable = TRUE,
       use.packages = FALSE,
       scalebox = 1.0)

#### table 5 effect of treatments on compliance forecasts ####

all <- all %>%
  mutate(
    pilot_prime = case_when(
      pilot_first == 1 & compliance_first == 1 ~ 1, # T1
      pilot_first == 0 & compliance_first == 1 ~ 0
      ), # T2
    compound = case_when(
      compliance_first == 0 ~ 1,
      compliance_first == 1 ~ 0
    )
  )

cphone_p <- lm_robust(ans_phone ~ pilot_prime, all)
cphone_c <- lm_robust(ans_phone ~ compound, all)

cquest_p <- lm_robust(ans_quest ~ pilot_prime, all)
cquest_c <- lm_robust(ans_quest ~ compound, all)

texreg(l = list(cphone_c, cquest_c),
       include.ci = FALSE,
       digits = 3,
       booktabs = TRUE,
       label="tab:compl-results",
       stars = c(0.01, 0.05, 0.1),
       caption = "Experimental effects for compliance forecasts",
       caption.above = TRUE,
       custom.header = list("\\textit{Dependent variable}:" = 1:2),
       custom.model.names = c("Compliance rate for call", "Compliance rate for question"),
       custom.coef.map = list("(Intercept)" = "Constant", "compound" = "Pilot results effect"),
       reorder.coef = c(2, 1),
       include.rsquared = F, include.adjrs = F, include.rmse = F, include.nobs = T,
       custom.note = "\n\\item %stars. OLS estimates with HC2 standard errors. \\textit{Compliance rate for call} refers to the percentage of treated RCT subjects who answer the MPA's call. \\textit{Compliance rate for question} refers to the percentage of treated RCT subjects answering the call who also answer the MPA's question.",
       float.pos = "h",
       threeparttable = TRUE,
       use.packages = FALSE,
       scalebox = 1.1)


#### table 6 ATEs conditional on familiarity with IVR technology ####

## regression models ##

fam_phone_c <- lm_robust(ans_phone ~ compound*familiar, data = all)
fam_phone_p <- lm_robust(ans_phone ~ pilot_prime*familiar, data = all)

fam_quest_c <- lm_robust(ans_quest ~ compound*familiar, data = all)
fam_quest_p <- lm_robust(ans_quest ~ pilot_prime*familiar, data = all)

texreg(l = list(fam_phone_c, fam_quest_c),
       include.ci = FALSE,
       digits = 3,
       booktabs = TRUE,
       label="tab:interact",
       stars = c(0.01, 0.05, 0.1),
       caption = "Effect heterogeneity for compliance forecasts",
       caption.above = TRUE,
       custom.header = list("\\textit{Dependent variable}:" = 1:2),
       custom.model.names = c("Compliance rate for call", "Compliance rate for question"),
       custom.coef.map = list("(Intercept)" = "Constant", "compound" = "Pilot results",
                              "familiar" = "Familiar",
                              "compound:familiar" = "Pilot results x Familiar"),
       reorder.coef = c(4, 2, 3, 1),
       include.rsquared = F, include.adjrs = F, include.rmse = F, include.nobs = T,
       custom.note = "\n\\item %stars. OLS estimates with HC2 standard errors. \\textit{Compliance rate for call} refers to the percentage of treated RCT subjects who answer the MPA's call. \\textit{Compliance rate for question} refers to the percentage of treated RCT subjects answering the call who also answer the MPA's question. \\textit{Familiar} refers to subjects who score themselves 'familiar' or 'very familiar' with ICT in governance.",
       float.pos = "h",
       threeparttable = TRUE,
       use.packages = FALSE,
       scalebox = 1.0)

#### figure 2 interaction plot ####

cf1 <- lm_robust(ans_quest ~ compound, data = all, subset = familiar == 1)
cf0 <- lm_robust(ans_quest ~ compound, data = all, subset = familiar == 0)

intmodel <- lm(ans_quest ~ compound*familiar, data = all)

intplot1 <- interact_plot(intmodel, pred = compound, modx = familiar, plot.points = TRUE,
              jitter = 0.1, point.shape = TRUE)

#ggsave("intplot1.pdf", intplot1, width = 8, height = 6)

simslopes <- sim_slopes(intmodel, pred = compound, modx = familiar, johnson_neyman = FALSE,
                        robust = "HC2", cond.int = TRUE)
#plot(simslopes)

## compare group means incl CIs
emm <- emmeans(intmodel, ~ compound | familiar)
intdf <- as_tibble(plot(emm, plotit = FALSE))

intdf <- intdf %>%
  mutate(compound = as_factor(compound),
         familiar2 = ifelse(familiar == 1, "Familiar", "Unfamiliar")
  )

data.segm1 <- data.frame(x = "0", y = intdf$the.emmean[intdf$compound == 0 & intdf$familiar == 1], xend = "1",
                         yend = intdf$the.emmean[intdf$compound == 1 & intdf$familiar == 1],
                         familiar2 = "Familiar", compound = "1")

data.segm2 <- data.frame(x = "0", y = intdf$the.emmean[intdf$compound == 0 & intdf$familiar == 0], xend = "1",
                         yend = intdf$the.emmean[intdf$compound == 1 & intdf$familiar == 0],
                         familiar2 = "Unfamiliar", compound = "1")

text1 <- data.frame(x = "1", y = 0.595, label = "Diff = -0.004 (p = 0.93)",
                    familiar2 = "Familiar", compound = "1")

text2 <- data.frame(x = "1", y = 0.595, label = "Diff = 0.16 (p < 0.001)",
                    familiar2 = "Unfamiliar", compound = "1")

plot <- ggplot(data = intdf, mapping = aes(x = compound, y = the.emmean)) +
  geom_point() +
  geom_errorbar(mapping = aes(x = compound, ymin = lower.CL, ymax = upper.CL), width = 0.1) +
  xlab("Compound information treatment") +
  ylab("Mean compliance forecast") +
  facet_wrap(~familiar2) +
  theme_pubclean()

intplot <- plot + geom_segment(data = data.segm1,
                    aes(x = x, y = y, yend = yend, xend = xend), inherit.aes = FALSE,
                    linetype = "dashed") +
  geom_segment(data = data.segm2,
               aes(x = x, y = y, yend = yend, xend = xend), inherit.aes = FALSE,
               linetype = "dashed") +
  geom_text(data = text1,
            aes(x = x, y = y, label = label), inherit.aes = FALSE) +
  geom_text(data = text2,
            aes(x = x, y = y, label = label), inherit.aes = FALSE)

pmatrix <- pwpm(emm, flip = FALSE) # p-values

ggsave("intplot.pdf", intplot, width = 10, height = 8)


#### figure 3 RCT count analysis ####

df <- plyr::ldply(list_all, data.frame)

df <- df %>%
  mutate(
    journal2 = factor(case_when(
      journal == "ajps" ~ "AJPS",
      journal == "apsr" ~ "APSR",
      journal == "jop" ~ "JOP"), levels = c("APSR", "AJPS", "JOP")),
    result2 = ifelse(result == 0, "Null effect",
                     ifelse(result == 1, "Effect", NA))
  )


df %>%
  group_by(journal) %>%
  summarize(n = n())

plot1 <- df %>%
  group_by(journal2, rct) %>%
  summarize(n = n()) %>%
  mutate(prop = n / sum(n))

p1 <- ggplot(data = plot1 %>% filter(rct == 1),
             mapping = aes(x = journal2, y = prop)) +
  geom_bar(stat = "identity", width = 0.6, color = "blue", fill = rgb(0.1,0.4,0.5,0.7)) +
  xlab("") + ylab("Proportion of studies with RCTs out of all articles") +
  geom_text(aes(label = n), position = position_dodge(width = 0.9), vjust = -0.25) +
  theme_light(base_size = 14) +
  ggtitle("Proportion of RCTs")

plot2 <- df %>%
  filter(rct == 1) %>%
  group_by(journal2, result2) %>%
  summarize(n = n()) %>%
  mutate(
    prop = n / sum(n),
    result2 = factor(result2, levels = c("Effect", "Null effect"))
  )

p2 <- ggplot(data = plot2, mapping = aes(x = journal2, y = prop, fill = result2)) +
  geom_bar(stat = "identity", color = "blue", width = 0.6) +
  scale_fill_manual(values = c(rgb(0.1,0.4,0.5,0.7), rgb(0.2,0.3,0.6,0.8))) +
  xlab("") + ylab("Proportion of RCTs with(out) null results out of all RCTs") +
  geom_text(aes(label = n), position = position_stack(vjust = 0.95)) +
  theme_light(base_size = 14) +
  guides(fill = guide_legend(title = NULL)) +
  ggtitle("Proportion of RCT null results")

journals <- ggarrange(p1, p2, nrow = 1, ncol = 2, align = "h")

ggsave("rct_plot.pdf", journals, width = 10, height = 6)

#### table 7 expert status by position ####
CrossTable(all$familiar, all$position, prop.r=TRUE, prop.c=FALSE,
           prop.t=FALSE, prop.chisq=FALSE)

#### table 8 expert status by wave ####
CrossTable(all$familiar, all$wave, prop.r=TRUE, prop.c=FALSE,
           prop.t=FALSE, prop.chisq=FALSE)

#---- 2. APPENDIX ANALYSES ----

#### table a.1 effect of pilot prime on compliance forecasts ####
texreg(l = list(cphone_p, cquest_p),
       include.ci = FALSE,
       digits = 3,
       booktabs = TRUE,
       label="tab:compl-results-pilot",
       stars = c(0.01, 0.05, 0.1),
       caption = "Effects of pilot prime on compliance forecasts",
       caption.above = TRUE,
       custom.header = list("\\textit{Dependent variable}:" = 1:2),
       custom.model.names = c("Compliance rate for call", "Compliance rate for question"),
       custom.coef.map = list("(Intercept)" = "Constant", "pilot_prime" = "Pilot prime"),
       reorder.coef = c(2, 1),
       include.rsquared = F, include.adjrs = F, include.rmse = F, include.nobs = T,
       custom.note = "\n\\item %stars. OLS estimates from regressions of compliance forecasts on the pilot prime indicator with HC2 standard errors. \\textit{Compliance rate for call} refers to the percentage of treated RCT subjects who answer the MPA's call. \\textit{Compliance rate for question} refers to the percentage of treated RCT subjects answering the call who also answer the MPA's question.",
       float.pos = "h",
       threeparttable = TRUE,
       use.packages = FALSE,
       scalebox = 1.1)



#### table a.2 heterogeneity by optimism ####
all <- all %>%
  mutate(optimism_group = factor(optimism_group, levels = c("Pessimistic", "Neutral", "Optimistic")))

opt_phone_c <- lm_robust(ans_phone ~ compound*optimism_group, data = all)
opt_phone_p <- lm_robust(ans_phone ~ pilot_prime*optimism_group, data = all)

opt_quest_c <- lm_robust(ans_quest ~ compound*optimism_group, data = all)
opt_quest_p <- lm_robust(ans_quest ~ pilot_prime*optimism_group, data = all)

texreg(l = list(opt_phone_c, opt_quest_c),
       include.ci = FALSE,
       digits = 3,
       booktabs = TRUE,
       label="tab:interact-opt",
       stars = c(0.01, 0.05, 0.1),
       caption = "Effect heterogeneity compliance forecasts",
       caption.above = TRUE,
       custom.header = list("\\textit{Dependent variable}:" = 1:2),
       custom.model.names = c("Compliance rate for call", "Compliance rate for question"),
       custom.coef.map = list("(Intercept)" = "Constant", "compound" = "Pilot results",
                              "optimism_groupNeutral" = "Neutral", "optimism_groupOptimistic" = "Optimistic",
                              "compound:optimism_groupNeutral" = "Pilot results x Neutral",
                              "compound:optimism_groupOptimistic" = "Pilot results x Optimistic"),
       reorder.coef = c(5, 6, 4, 3, 2, 1),
       include.rsquared = F, include.adjrs = F, include.rmse = F, include.nobs = T,
       custom.note = "\n\\item %stars. OLS estimates with HC2 standard errors. \\textit{Compliance rate for call} refers to the percentage of treated RCT subjects who answer the MPA's call. \\textit{Compliance rate for question} refers to the percentage of treated RCT subjects answering the call who also answer the MPA's question. \\textit{Optimistic} refers to subjects who ranked themselves Very optimistic or Optimistic on a 1-5 scale that asked 'How optimistic are you about the potential for information technology to improve governance in Pakistan?'",
       float.pos = "h",
       threeparttable = TRUE,
       use.packages = FALSE,
       scalebox = 1.0)


#### table a.3 effect of treatments on forecasting error for compliance ####
all <- all %>%
  mutate(compl_error_phone = abs(0.731 - ans_phone),
         compl_error_quest = abs(0.238 - ans_quest))

ep_pilot <- lm_robust(compl_error_phone ~ pilot_prime, data = all)
ep_compl <- lm_robust(compl_error_phone ~ compliance_first, data = all)

eq_pilot <- lm_robust(compl_error_quest ~ pilot_prime, data = all)
eq_compl <- lm_robust(compl_error_quest ~ compliance_first, data = all)

all$call_eval_mpa_dist <- abs(all$call_eval_mpa - 0.09) # initial error (0.09 is pilot estimate we told them about)

em_pilot <- lm_robust(call_eval_mpa_dist ~ pilot_first, data = all)
em_compl <- lm_robust(call_eval_mpa_dist ~ compliance_first, data = all)
em_cp <- lm_robust(call_eval_mpa_dist ~ compliance_first*pilot_first, data = all)

texreg(l = list(ep_pilot, ep_compl, eq_pilot, eq_compl, em_pilot, em_compl, em_cp),
       include.ci = FALSE,
       digits = 3,
       booktabs = TRUE,
       label="tab:error",
       stars = c(0.01, 0.05, 0.1),
       caption = "Effects on Forecasting Errors",
       caption.above = TRUE,
       custom.header = list("Phone" = 1:2, "Question" = 3:4, "Call MPA" = 5:7),
       custom.model.names = c("Compl.", "Compl.",
                              "Compl.", "Compl.",
                              "ITT", "ITT", "ITT"),
       custom.coef.map = list("pilot_prime" = "Pilot prime", "compliance_first" = "Compliance prime","compliance_first:pilot_first" = "Interaction PxC", "pilot_first" = "Pilot prime", "(Intercept)" = "Constant"),
       include.rsquared = F, include.adjrs = F, include.rmse = F, include.nobs = T,
       custom.note = "\n\\item %stars. OLS estimates from regressions of compliance and ITT forecast errors on treatment group indicators with HC2 standard errors. Errors are defined as the absolute distance of initial forecasts vs actual compliance rates (73.1\\% for answering phone, 23.8\\% for answering question conditional on picking up the phone) or vs pilot information we told subjects about in the treatment (0.09 sd).",
       float.pos = "h",
       threeparttable = TRUE,
       use.packages = FALSE,
       scalebox = 0.9)


#### table a.4 effects of treatments on updating direction ####

## EVAL MPA 
all$call_eval_mpa_update <- all$call_eval_mpa_rev - all$call_eval_mpa # revised - initial estimate
all$call_eval_mpa_update_right_dir <- all$call_eval_mpa_update * ifelse(all$call_eval_mpa < 0.09, 1, -1)
## estimate is unit change toward the right direction

um_pilot <- lm_robust(call_eval_mpa_update_right_dir ~ pilot_first, data = all)
um_compl <- lm_robust(call_eval_mpa_update_right_dir ~ compliance_first, data = all)
um_cp <- lm_robust(call_eval_mpa_update_right_dir ~ compliance_first*pilot_first, data = all)

## EVAL GOV
all$call_eval_gov_update <- all$call_eval_gov_rev - all$call_eval_gov # revised - initial estimate
all$call_eval_gov_update_right_dir <- all$call_eval_gov_update * ifelse(all$call_eval_gov > -0.05, 1, -1)

ug_pilot <- lm_robust(call_eval_gov_update_right_dir ~ pilot_first, data = all)
ug_compl <- lm_robust(call_eval_gov_update_right_dir ~ compliance_first, data = all)
ug_cp <- lm_robust(call_eval_gov_update_right_dir ~ compliance_first*pilot_first, data = all)

## ACCOUNTABILITY
all$call_account_update <- all$call_account_rev - all$call_account # revised - initial estimate
all$call_account_update_right_dir <- all$call_account_update * ifelse(all$call_account < 0.12, 1, -1)

ua_pilot <- lm_robust(call_account_update_right_dir ~ pilot_first, data = all)
ua_compl <- lm_robust(call_account_update_right_dir ~ compliance_first, data = all)
ua_cp <- lm_robust(call_account_update_right_dir ~ compliance_first*pilot_first, data = all)


texreg(l = list(um_pilot, um_compl, um_cp, ug_pilot, ug_compl, ug_cp, ua_pilot, ua_compl, ua_cp),
       include.ci = FALSE,
       digits = 3,
       booktabs = TRUE,
       label="tab:update-dir",
       stars = c(0.01, 0.05, 0.1),
       caption = "Effects on Correct Updating Direction",
       caption.above = TRUE,
       custom.header = list("MPA" = 1:3, "Government" = 4:6, "Accountability" = 7:9),
       custom.model.names = c("ITT", "ITT", "ITT",
                              "ITT", "ITT", "ITT",
                              "ITT", "ITT", "ITT"),
       custom.coef.map = list("pilot_first" = "Pilot prime", "compliance_first" = "Compliance prime", "compliance_first:pilot_first" = "Interaction PxC", "(Intercept)" = "Constant"),
       include.rsquared = F, include.adjrs = F, include.rmse = F, include.nobs = T,
       custom.note = "\n\\item %stars. OLS estimates from regressions of a continuous variable measuring the extent of updating in the correct direction on treatment group indicators with HC2 standard errors. Updating in the correct direction measures the unit change in the update toward the right estimate for each outcome. Depending on where the subjects started (above or below the pilot estimate they were informed about), positive values indicate moving to the correct direction from \\textit{either} side.",
       float.pos = "h",
       threeparttable = TRUE,
       use.packages = FALSE,
       scalebox = 0.8)


#### table a.5 effect of treatments on update magnitude ####

## EVAL MPA
magmpa_pilot <- lm_robust(call_eval_mpa_update ~ pilot_first, data = all)
magmpa_compl <- lm_robust(call_eval_mpa_update ~ compliance_first, data = all)
magmpa_cp <- lm_robust(call_eval_mpa_update ~ compliance_first*pilot_first, data = all)


## EVAL GOV
maggov_pilot <- lm_robust(call_eval_gov_update ~ pilot_first, data = all)
maggov_compl <- lm_robust(call_eval_gov_update ~ compliance_first, data = all)
maggov_cp <- lm_robust(call_eval_gov_update ~ compliance_first*pilot_first, data = all)

## ACCOUNTABILITY
magac_pilot <- lm_robust(call_account_update ~ pilot_first, data = all)
magac_compl <- lm_robust(call_account_update ~ compliance_first, data = all)
magac_cp <- lm_robust(call_account_update ~ compliance_first*pilot_first, data = all)

texreg(l = list(magmpa_pilot, magmpa_compl, magmpa_cp, maggov_pilot, maggov_compl, maggov_cp, magac_pilot, magac_compl, magac_cp),
       include.ci = FALSE,
       digits = 3,
       booktabs = TRUE,
       label="tab:update-mag",
       stars = c(0.01, 0.05, 0.1),
       caption = "Effects on Updating Magnitude",
       caption.above = TRUE,
       custom.header = list("MPA" = 1:3, "Government" = 4:6, "Accountability" = 7:9),
       custom.model.names = c("ITT", "ITT", "ITT",
                              "ITT", "ITT", "ITT",
                              "ITT", "ITT", "ITT"),
       custom.coef.map = list("pilot_first" = "Pilot prime", "compliance_first" = "Compliance prime", "compliance_first:pilot_first" = "Interaction PxC", "(Intercept)" = "Constant"),
       include.rsquared = F, include.adjrs = F, include.rmse = F, include.nobs = T,
       custom.note = "\n\\item %stars. OLS estimates from regressions of the extent of updating of ITT forecasts on treatment group indicators with HC2 standard errors. Updating is defined as the difference between revised and initial ITT estimate, respectively for all outcomes.",
       float.pos = "h",
       threeparttable = TRUE,
       use.packages = FALSE,
       scalebox = 0.7)


#### table a.6 effect of treatments on forecasts - call only ####

texreg(l = list(call_mpa_p, call_mpa_c, call_mpa_pc,
                 call_gov_p, call_gov_c, call_gov_pc,
                 call_acc_p, call_acc_c, call_acc_pc),
       include.ci = FALSE,
       digits = 3,
       booktabs = TRUE,
       label="tab:full-call",
       stars = c(0.01, 0.05, 0.1),
       caption = "Effects on Forecasted ITT (Call only treatment)",
       caption.above = TRUE,
       custom.header = list("MPA" = 1:3, "Government" = 4:6, "Accountability" = 7:9),
       custom.model.names = c("ITT", "ITT", "ITT",
                              "ITT", "ITT", "ITT",
                              "ITT", "ITT", "ITT"),
       custom.coef.map = list("pilot_first" = "Pilot prime", "compliance_first" = "Compliance prime", "pilot_first:compliance_first" = "Interaction PxC", "(Intercept)" = "Constant"),
       include.rsquared = F, include.adjrs = F, include.rmse = F, include.nobs = T,
       custom.note = "\n\\item %stars. OLS estimates from regressions of all forecasted outcomes on treatment group indicators with HC2 standard errors. We focus here on the forecasts for the effects of the Call only (generic) treatments on the individual indices separately, i.e. forecasted ITT for MPA evaluation, Government evaluation and Prospects for accountability.",
       float.pos = "h",
       threeparttable = TRUE,
       use.packages = FALSE,
       scalebox = 0.75)

#### table a.7 effect of treatments on forecasts - responsive only ####

texreg(l = list(resp_mpa_p, resp_mpa_c, resp_mpa_pc,
                 resp_gov_p, resp_gov_c, resp_gov_pc,
                 resp_acc_p, resp_acc_c, resp_acc_pc),
       include.ci = FALSE,
       digits = 3,
       booktabs = TRUE,
       stars = c(0.01, 0.05, 0.1),
       label="tab:full-resp",
       caption = "Effects on Forecasted ITT (Responsive treatment)",
       caption.above = TRUE,
       custom.header = list("MPA" = 1:3, "Government" = 4:6, "Accountability" = 7:9),
       custom.model.names = c("ITT", "ITT", "ITT",
                              "ITT", "ITT", "ITT",
                              "ITT", "ITT", "ITT"),
       custom.coef.map = list("pilot_first" = "Pilot prime", "compliance_first" = "Compliance prime", "pilot_first:compliance_first" = "Interaction PxC", "(Intercept)" = "Constant"),
       include.rsquared = F, include.adjrs = F, include.rmse = F, include.nobs = T,
       custom.note = "\n\\item %stars. OLS estimates from regressions of all forecasted outcomes on treatment group indicators with HC2 standard errors. We focus here on the forecasts for the effects of the Responsive treatment on the individual indices separately, i.e. forecasted ITT for MPA evaluation, Government evaluation and Prospects for accountability.",
       float.pos = "h",
       threeparttable = TRUE,
       use.packages = FALSE,
       scalebox = 0.75)


#### table c.1 summary statistics ####
all$familiarity <- as.numeric(all$familiarity)

all2 <- all %>%
  select(familiarity, optimism, call_avg, call_account, call_eval_gov, call_eval_mpa,
         resp_avg, responsive_account, responsive_eval_gov, responsive_eval_mpa,
         ans_phone, ans_quest)

col_summary <- function(df, fun) {
  out <- vector("double", length(df))
  for (i in seq_along(df)) {
    out[i] <- fun(df[[i]], na.rm = TRUE)
  }
  out
}

sum_stats <- tibble(
  Variable = colnames(all2),
  Mean = round(col_summary(all2, mean), 3),
  `St. Dev.` = round(col_summary(all2, sd), 3),
  Min = round(col_summary(all2, min), 2),
  Max = round(col_summary(all2, max), 2)
) %>%
  mutate(Variable = c("Familiarity", "Optimism", "Forecast Call Avg.", "Forecast Call Accountability", "Forecast Call Government", "Forecast Call MPA", "Forecast Responsive Avg.", "Forecast Responsive Accountability", "Forecast Responsive Government", "Forecast Responsive MPA", "Compliance rate phone", "Compliance rate question"))

stargazer(sum_stats, type = "latex", summary = FALSE, rownames = FALSE,
          title = "Summary statistics on main variables",
          header = FALSE,
          label="tab:sum-stats"
)

#### table c.2 experimental subjects ####

respondents <- all %>%
  group_by(wave) %>%
  summarize(n = n()) %>%
  bind_rows(tibble(wave = "total", n = 280))

stargazer(respondents, type = "latex", summary = FALSE, rownames = FALSE,
          title = "Experimental Subjects",
          header = FALSE,
          label="tab:sample")

#### table d.1 balance test for familiarity and optimism ####

all <- all %>%
  filter(compliance_first == 1 | compliance_first == 0) %>%
  mutate(
    pilot_prime = case_when(
      pilot_first == 1 & compliance_first == 1 ~ 1,
      pilot_first == 0 & compliance_first == 1 ~ 0
    ),
    t1 = case_when(
      compliance_first == 1 & pilot_first == 1 ~ 1,
      TRUE ~ 0
    ),
    t2 = case_when(
      compliance_first == 1 & pilot_first == 0 ~ 1,
      TRUE ~ 0
    ),
    t3 = case_when(
      compliance_first == 0 & pilot_first == 1 ~ 1,
      TRUE ~ 0
    ),
    t4 = case_when(
      compliance_first == 0 & pilot_first == 0 ~ 1,
      TRUE ~ 0
    )
  )

fam_ftest <- lm_robust(familiarity ~ t2 + t3 + t4, all) #t1 omitted dummy
fam_bal1 <- lm_robust(familiarity ~ compliance_first, all)
fam_bal2 <- lm_robust(familiarity ~ pilot_first, all)

opt_ftest <- lm_robust(optimism ~ t2 + t3 + t4, all) # t1 omitted dummy
opt_bal1 <- lm_robust(optimism ~ compliance_first, all)
opt_bal2 <- lm_robust(optimism ~ pilot_first, all)

texreg(l = list(fam_ftest, fam_bal1, fam_bal2,
                opt_ftest, opt_bal1, opt_bal2),
       include.ci = FALSE,
       digits = 3,
       booktabs = TRUE,
       label="tab:bal1",
       stars = c(0.01, 0.05, 0.1),
       caption = "Balance tests",
       caption.above = TRUE,
       custom.header = list("Familiarity" = 1:3, "Optimism" = 4:6),
       custom.model.names = c("(1)", "(2)", "(3)", "(4)", "(5)", "(6)"),
       custom.coef.map = list("(Intercept)" = "Mean T1", "t2" = "T2", "t3" = "T3", "t4" = "T4", "compliance_first" = "Compliance prime", "pilot_first" = "Pilot prime"),
       include.rsquared = F, include.adjrs = F, include.rmse = F,
       include.nobs = T, include.fstatistic = T,
       custom.note = "\n\\item %stars. OLS estimates from regressions of familiarity [1, 4] and optimism [1, 5] variables on treatment group indicators with HC2 standard errors. For Compliance first (T1, T2), the reference group is (T3, T4) and for Pilot prime (T1, T3), the reference group is (T2, T4).",
       float.pos = "h",
       threeparttable = TRUE,
       use.packages = FALSE,
       scalebox = 1.0)

#### table d.2 balance test I for forecaster position ####

all$position <- as_factor(all$position)
pos_bal1 <- multinom(position ~ compliance_first, all,
                     trace = FALSE)
pos_c <- tidy(pos_bal1)

pos_bal2 <- multinom(position ~ pilot_first, all,
                     trace = FALSE)
pos_p <- tidy(pos_bal2)

pos_global <- multinom(position ~ t2 + t3 + t4, all,
                       trace = FALSE)
pos_t <- tidy(pos_global)

stargazer(pos_bal1,
          type = "latex",
          label="tab:bal2",
          #align = TRUE,
          dep.var.labels=c("Faculty","Other research staff", "Postdoc"),
          covariate.labels=c("Compliance prime (T1+T2)", "Mean T3+T4"),
          #keep.stat = "n",
          no.space = TRUE,
          title = "Balance Test I For Forecaster Position",
          header = FALSE,
          notes.append = TRUE, notes = "\\parbox[t]{10cm}{Multinomial logit coefficients from regressions of forecaster position with omitted baseline category Graduate students.}",
          notes.align = "r"
)

#### table d.3 balance test II for forecaster position ####

stargazer(pos_bal2,
          type = "latex",
          label="tab:bal3",
          #align = TRUE,
          dep.var.labels=c("Faculty","Other research staff", "Postdoc"),
          covariate.labels=c("Pilot prime (T1+T3)", "Mean T2+T4"),
          #keep.stat = "n",
          no.space = TRUE,
          title = "Balance Test II For Forecaster Position",
          header = FALSE,
          notes.append = TRUE, notes = "\\parbox[t]{10cm}{Multinomial logit coefficients from regressions of forecaster position with omitted baseline category Graduate students.}",
          notes.align = "r"
)

#### table d.4 balance test III for forecaster position ####

stargazer(pos_global,
          type = "latex",
          label="tab:bal4",
          #align = TRUE,
          dep.var.labels=c("Faculty","Other research staff", "Postdoc"),
          covariate.labels=c("T1", "T2", "T3", "Mean T1"),
          #keep.stat = "n",
          no.space = TRUE,
          title = "Balance Test III For Forecaster Position",
          header = FALSE,
          notes.append = TRUE, notes = "\\parbox[t]{10cm}{Multinomial logit coefficients from regressions of forecaster position with omitted baseline category Graduate students.}",
          notes.align = "r"
)



